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Recently, we have derived an effective Cahn-Hilliard equation for the phase separation dynamics of 
active Brownian particles by performing a weakly non-linear analysis of the effective hydrodynamic 
equations for density and polarization [Phys. Rev. Lett. 112, 218304 (2014)]. Here we develop 
and explore this strategy in more detail and show explicitly how to get to such a large-scale, mean- 
field description starting from the microscopic dynamics. The effective free energy emerging from 
this approach has the form of a conventional Ginzburg-Landau function. On the coarsest scale, 
our results thus agree with the mapping of active phase separation onto that of passive fluids with 
attractive interactions through a global effective free energy (mobility-induced phase transition). 

Particular attention is paid to the square-gradient term necessary for the dynamics. We finally 
discuss results from numerical simulations corroborating the analytical results. 
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I. INTRODUCTION 

The separation of a suspension of passive colloidal par¬ 
ticles into liquid and vapor is a complex and rather well- 
studied phenomenon [T|. For phase separation to occur, 
sufficiently strong attractions between particles have to 
be present. This is well understood from the perspective 
of thermodynamics: The potential energy gained by the 
suspension forming the dense phase compensates for the 
loss of entropy. 

Recently, a separation into dense and dilute regions has 
also been reported for “active” suspensions composed of 
self-propelled colloidal particles HHU- For a brief per¬ 
spective on this phenomenon see Ref. [5] and Refs. H5H2] for 
more general recent reviews on various aspects of active 
matter. By constantly converting external energy into 
directed motion, such systems of self-propelled particles 
can be driven into a non-equilibrium steady state. Steady 
dynamic states of orientationally ordered collective mo¬ 
tion can arise [TOHIOj , but also swirling and turbulent-like 
situations were observed [20M25] . There are different pos¬ 
sibilities how to provide the energy input in experiments: 
light, if sufficiently strong, can create a temperature gra¬ 
dient leading to self-thermophoresis [25] ;25j, or lead to 
the local demixing of a binary solvent |26| . In most ex¬ 
periments, however, energy is provided chemically, e.g., 
through the decomposition of hydrogen peroxide [27: or 
the release of stored ions [25]. All of these mechanisms 
are based on a symmetry breaking on the particle level. 
In the case of spherical Janus particles [25 it is provided 
by different surface properties of typically two distinct 
hemispheres. 

What is intriguing from a fundamental perspective is 
that the phase separation of active colloids strongly re¬ 
sembles the phase separation in a passive suspension, but 
occurs also for purely repulsive particles. This has been 
demonstrated convincingly in computer simulations of 
a minimal model by a number of groups [30)433j . This 


minimal model of active Brownian particles incorporates 
the two basic physical ingredients: volume exclusion and 
persistence of motion, i.e., particles interact via repul¬ 
sive potentials (or hard-core exclusion) and every parti¬ 
cle has an orientation along which it “swims” at constant 
speed. The particle orientations evolve independently 
and stochastically. The microscopic reason for particle 
accumulation is simple (see, e.g., the kinetic model in 
Ref. ITTT[) : Due to the persistence of the self-propelled mo¬ 
tion and the excluded volume, particles block each other 
on the time scale that it takes for orientations to decorre¬ 
late. If the mean collision time is shorter (the suspension 
is sufficiently dense) clustering ensues. Simulations with 
idealized boundary conditions show that hydrodynamic 
interactions due to the solvent can modify these time 
scales [551 155] . On the other hand, experiments |T| in¬ 
dicate that for colloidal swimmers the basic scenario is 
robust, and we will neglect hydrodynamic interactions in 
the following. 

Even if the microscopic mechanism is known, the col¬ 
lective large-scale and phase behavior is still highly non¬ 
trivial. Tailleur and Cates have been the first to realize 
that the phenomenon of phase separation in active sys¬ 
tems can be explained by the dependence of the parti¬ 
cle mobility on local density [5H1 - !55| . The resulting the¬ 
ory is referred to as “mobility induced phase separation” 
(MIPS). Within this framework it has been demonstrated 
that the temporal evolution of the coarse-grained den¬ 
sity can be mapped onto that of an effective bulk free 
energy. Good agreement between particle-resolved sim¬ 
ulations of active Brownian particles and the numerical 
solution of the coarse-grained density has been demon¬ 
strated jUlED]- In analogy to the mean-field treatment 
of phase separation in passive suspensions (with attrac¬ 
tive forces), binodal and spinodal lines are identified from 
the minima and inflection points of the effective bulk free 
energy, respectively, from which a schematic phase dia¬ 
gram is constructed. Both lines merge at a single critical 
point. The resulting free energy involves only the density, 
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a result that is based on the elimination of the polariza¬ 
tion through neglecting temporal and spa tial derivatives 
(a detailed discussion follows in Sec. III). In particular, 
this treatment produces a bulk free energy that does not 
contain a term that penalizes sharp interfaces. It is thus 
not able to describe the dynamics of domain formation 
and coarsening. In order to cure this shortcoming, Cates 
and coworkers have argued that the dynamics has to be 
augmented by gradient terms that are not derivable from 
a free energy 


The purpose of the present paper is to follow the com¬ 
plete path from the microscopic dynamics to the large- 
scale Cahn-Hilliard equation [32] and to give a compre¬ 
hensive derivation of the results obtained in Ref. 03] To 
this end we derive effective hydrodynamic equations and 
perform a weakly non-linear analysis |44j . Our results 
refine the MIPS scenario in two ways: (i) Our system¬ 
atic derivation leads to a square-gradient term that com¬ 
pletes the bulk free energy to the conventional Ginzburg- 
Landau form, (ii) We clarify the validity of a description 
in terms of the density alone and show that it can strictly 
hold only close to the loss of linear stability. The physi¬ 
cal reason is that this instability occurs on length scales 
that are larger than the persistence length of the directed 
motion and thus allow us to study the coarse-grained den¬ 
sity alone. The microscopic “trapping” of particles due 
to their directed motion then enters as an effective at¬ 
traction. Quenching deeper into the two-phase region, 
the non-equilibrium nature of active suspensions will be¬ 
come evident again and the coupled evolution of density 
and orientation has to be considered. 


Our paper is organized as follows: In Sec. [TTj, we dis¬ 
cuss active Brownian particles as a minimal model for 
self-propelled disks. Starting from the full many-body 
dynamics, we sketch the derivation of the effective hy¬ 
drodynamic equations of motion following Ref. 1321 In 
Sec. |III[ we explore the consequences of the adiabatic ap¬ 
proximation of the hydrodynamic equations leading to 
an effective equilibrium theory for the density alone. We 
then describe the weakly non-linear analysis in Sec. |IV[ 
the ramifications of which are discussed in Sec. El Con¬ 
clusions and outlook are given in Sec. m 


II. MEAN-FIELD THEORY 


A. Minimal model for self-propelled disks 

The model we study consists of N identical, interact¬ 
ing disks with free diffusion coefficient Dq moving in a 
periodic box with area A. Each disk has an orientation 
e*, = (cos <pk, sin (pk) T that undergoes free rotational dif¬ 
fusion with diffusion coefficient D r . We consider both 
diffusion coefficients to be constant and eliminate them 
through rescaling time t K > t/D r and length r i—»• £r, 
where £ = yjD q / D r . Neglecting hydrodynamic interac¬ 


tions, the coupled equations of motion become 

r fc = -X7U + v 0 e k +£ fc , (1) 

where U = J2k<k' u ( l r fc — r k'\) the potential energy 
(in units of the thermal energy) with pair potential zt(r), 
and the noise correlations due to the solvent read 

(t k m T k'(t')) = 2l6 kkl 6(t-t'). (2) 

The orientational angles obey 

(■')) = 2 5 k k'S(t - t'). (3) 

We are thus left with two dimensionless parameters defin¬ 
ing a non-equilibrium state point: the number density 
p = Nl 2 /A and the free speed vq- 


B. Derivation 


An equivalent representation of the equations of mo¬ 
tion Eq. 0 is given by the Smoluchowski equation 

N ( Q2 

dtV’iv = ^2 | V*; • [(Vfc/7) - v 0 e k + Vfc] + j Vw 

(4) 

for the joint probability V’Ar({rfc, fk}', t) of particle posi¬ 
tions and their orientations. This representation is conve¬ 
nient since it allows for systematic approximations. Fol¬ 
lowing Ref. [32], we aim to derive a closed equation of 
motion for the one-point particle density 


Vh(iT,<£>i ;t) = 


dr 2 


• - dr at 



• • • difiN Ni/j n 


(5) 


integrating out all other particles. The local density 


P{r,t) 


dip V>i(r,^;t) 


( 6 ) 


corresponds to the probability of finding a particle at 
position r at time t. 

Inserting Eq. Q into the time derivative of Eq. ([5|, 
we obtain 

dt'ipi = -V • [F + u 0 e^i - VV’i] + d^ipi (7) 

dropping the particle index. The force F(r, ip;t) is due 
to interactions of a (fixed) particle with its surrounding 
particles averaged over their accessible positions and ori¬ 
entations. It thus couples to the two-point density 


ip 2 {ri,ipi,r 2 ,ip2;t) = J dr 3 ---drjv 

x J dip 3 -■ ■ d<p N N(N - (8) 


Even in a homogeneous suspension of active particles, the 
force F(r, ip\t) does not vanish. Rather, there is a force 
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FIG. 1: Sketch of the force decomposition. The total aver¬ 
age force F is approximately decomposed into a component 
along the particle orientation e and along the local one-point 
density gradient Vtpi. Also shown is one of the surrounding 
particles at distance r enclosing the angle 9 with the particle 
orientation. 


imbalance on each particle due to the directional motion, 
which implies a higher probability to find another particle 
in front of it (looking along its orientation) than behind. 
One finds F = e • F = —pCV'i with the force imbalance 
quantified through the coefficient 132 

n oo p2ir 

C = dr r[— u(r)} / d# cos 8g{r,0) (9) 

Jo Jo 

with pair distribution function g(r,9). Here, 8 is the 
angle between the orientation of a particle and the dis¬ 
placement vector to another particle at a distance r, see 
Fig. [l] For an inhomogeneous density p{ r) also £ de¬ 
pends in principle on the position r. Close to the dynam¬ 
ical instability of the homogeneous density profile we can 
neglect this spatial dependence and in the following we 
assume £ = £(p, Vq) to be a state function. 

We decompose the force 


F = (e ■ F)e + <5F ss (e ■ F)e + Fy Vifti 


( 10 ) 


into a component along the particle orientation due to 
the force imbalance and a component 5F perpendicular. 
The later describes an “evasive” motion leading to an 
effective diffusion. Hence, as a closure we only keep the 
projection of <5F onto the density gradient Vip i. This is a 
good approximation as long as | <5F | ~ 1 is much smaller 
compared to p( ~ v 0 1 [cf. Eq. (41)]. Rearranging 
V0i • F leads to the formal expression 1 "^ 


F, 


[V^i - (e ■ VVb)e] ■ F 

|VV’i | 2 


( 11 ) 


Inserting Eq. (101 into Eq. (J7|, the mean-field evolution 
equation for the joint probability of position and orien¬ 
tation thus reads 


d t ^i = -V • [v(p)e - D e V]V>i + (12) 


1 There is a typo in Ref. m below Eq. (11) (it should read | V"0i | 2 
in the expression for D). Consequently, Eq. H} is not an ex¬ 
pansion in the density gradient but follows from the different 
magnitudes, e • F . 


with effective diffusion coefficient D e = 1 — Fy >0 and 
effective speed 


v(p) = v 0 - pC 


(13) 


depending on the local density. 

For a homogeneous suspension, p{ r, t) = p is a con¬ 
stant. Assuming D e also to be constant, Eq. (121 for¬ 


mally corresponds to the stochastic evolution of a single 
self-propelled particle. We can then calculate the mean- 
square displacement and from that the long-time diffu¬ 
sion coefficient H51 


D(p,v 0 ) = D e (p) + ^[v(p)} 2 . 


(14) 


This relation has indeed been confirmed in several numer¬ 
ical studies of active suspensions [32:, [39], and therefore 
in the remainder we will treat D e (p) as spatially uni¬ 
form and independent of the speed. This constitutes our 
final approximation closing Eq. (12). Clearly, D e corre¬ 


sponds to the diffusion coefficient of the passive suspen¬ 
sion (^0 = 0). 


To briefly conclude, Eq. (12) describes the evolution 


of the active suspension on a coarse-grained level, into 
which the effects of microscopic particle interactions en¬ 
ter through two effective parameters: D e and the force 
imbalance £. Within the theory, every state point (p,v o) 
is fully characterized by these two parameters. However, 
the theory cannot make predictions about their values, 
for which we would have to make further assumptions or 
measure them in particle-resolved computer simulations 
(as has been done in Ref. [ 


C. Hydrodynamic equations 


For the evolution of the local density Eq. § one finds 
dtp = -V • Kp)p - D C V P \ (15) 


using Eq. (12). This equation expresses number conser¬ 


vation with a particle current that is given by a diffusive 
term —D c \7p and a current vp proportional to the polar¬ 
ization or orientational order parameter 


P(M) = J dl P (16) 

For p / 0, particles in a coarse-grained volume have a 
preferred orientation leading to a net particle current. 
This orientation evolves according to 

dtp = —'VP(p) + -D c V 2 p - P (17) 


with pressure 


p (p) = 2 V (P)P 


(18) 


resulting from the directed motion of the particles. In¬ 
serting Eq. (13), for sufficiently large £ this pressure 
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wave vector q 

FIG. 2: Dispersion relation a(q) for two speeds vo below and 
above the critical speed v c — 1.15 corresponding to = 1 
and v* = 1. For vo > v c large-scale perturbations ( i.e small 
wave vectors q > 0) become unstable. 


becomes a non-monotonic function of density. Hence, 
for local densities p > vq/{2Q the density gradient, 
VP(p) = P'(p)Vp, changes sign so that particles migrate 
towards denser regions. The second term in Eq. © is 
akin to a “viscosity term” and the last term describes the 
local relaxation due to the rotational diffusion. 

While we have derived Eqs. (15) and © starting 
from the full microscopic Smoluchowski equation Qj the 
formal structure of the result is not surprising since it 
has to reflect macroscopic conservation laws and symme¬ 
tries. Note that run-and-tumble dynamics leads to the 
same hydrodynamic equations (cf. Ref. [37) and see Ap¬ 
pendix |0. Although colloidal self-propelled particles of 
course move in a solvent, the decision to neglect hydro- 
dynamic interactions places the resulting effective hydro- 
dynamic theory into the field of what is sometimes called 
“dry active matter” [7|. Similar equations are, e.g., ob¬ 
tained in the Toner-Tu continuum treatment El of polar 
active systems [46] (47 ]. In that case the alignment of ori¬ 
entations leads to nonlinear terms of the polarization p 
in Eq. ([T7| and thus to dynamical collective behavior. 


D. Linear stability 


Clearly, any constant density with p = 0 is a station¬ 
ary solution of Eqs. dl5|) and ([T7|. To study the stability 


of the homogeneous density p with respect to small per¬ 
turbations, we set p = p + Sp and rewrite the equations 
of motion (15) and © 


d t 6p = ~a\7 -p + D e V 2 Sp + CV - (pip), (19) 
d t p = -/3\75p + P> c V 2 p - p + (SpVSp. (20) 


Here, we have separated the non-linear terms and defined 
the coefficients 

a = v{p) = v 0 - p(, P=^(v 0 -2pC). (21) 

Dropping the non-linear terms and inserting solutions of 
the form 5p, p oc e cr ^ 9 - |t+lq ' r , we obtain the dispersion 


relation 

<r(q) = -\- D eq 2 + \ V 1 - 4a/3g 2 

= (P e + a/3)q 2 - ( a/3) 2 q 4 + 0(q 6 ), 


( 22 ) 


which quantifies the growth rate of a perturbation with 
wave vector q (see Fig. 0. Hence, for D e + a(3 < 0 
a smooth perturbation of the homogeneous density on 
small q does not decay anymore but grows, leading to a 
dynamical instability. Solving this condition implies an 
instability line v c {p) of critical speeds (for simplicity we 
only consider the smaller solution) such that for vq > v c 
the suspension becomes linearly unstable. The growth 
rate <j(q) is maximized for the wave length 


2 _ 1 / 

Qo ~ 4 Dl) ’ 


(23) 


which will thus dominate the morphology at early stages 
after the onset of the instability. 


III. ADIABATIC APPROXIMATION 


For slowly varying fields, we neglect the temporal 
derivative as well as the viscosity term in Eq. ( [Tt| to 
obtain the adiabatic solution 


P ~ Pad = -VP = --V(^p). 


(24) 


Inserting Eq. (13) for the density-dependent effective 


speed v{p) there are two possibilities to rewrite this ex¬ 
pression, which we now discuss. 


A. Relaxation of density 


First, we push the gradient to the right, 
Pad = ~\{vo ~ 2pC)Vp. 


(25) 


We have assumed that ( is constant, which holds for the 
homogeneous density profile. Inserting p a d into Eq. (15) 
leads to the simple evolution equation for the density 

dtp = V • [D(p)Vp] (26) 

with collective diffusion coefficient 

V(p) = D e +^(v 0 -pC)(vo-2pQ. (27) 

The global density at which the diffusion coefficient 
switches its sign is given by the condition V(p) = 0. 
It signals the onset of a dynamical instability at which 
the homogeneous density profile p{ r, t) = p becomes (lin¬ 
early) unstable and small perturbations start to grow. 
Using the coefficients introduced in Eq. (21) we find 
V(p) = D e + a/3 and therefore cr(q) = —V(pjq 2 demon¬ 
strating that it is the same instability discussed in the 
previous section. 
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B. Free energy 

Alternatively, to illustrate the concept of an effective 
free energy as advocated by Cates and Tailleur ES1EZ1, 
we cast the evolution equation for the density into a form 
that involves the functional derivative of a potential func¬ 
tion, which thus can be interpreted as an effective free 
energy. This is achieved by pulling out the Nabla opera¬ 
tor, 


^ t'Pad = 2 V 


+ §<V 


(28) 


Inserting this result into Eq. (15), we now find the evo¬ 
lution equation 


d t p = -V • (up ad - D c Vp) = V 

implying the functional 

F ad [p\ = J d 2 r /ad(p(r)) 


2 $F ad 


5p 


(29) 


(30) 




with bulk free energy density 

/ad (p) = \( K D e + y) P 2 - yoCp 3 + yCV- (31) 

We will refer to this function as the adiabatic free energy 
density. It has the typical form of a Landau function 
often encountered in the study of critical phenomena and 
phase transitions )48j . 

In order to discuss the phase diagram following from 
the adiabatic solution, we rewrite the bulk free energy 
density as a symmetric function plus a linear term, 


FIG. 3: Schematic phase diagram for the adiabatic mean- 
field free energy density Eq. ( [32] ) using </(do) = (3Do)/(4p») 
with critical point at p„ = 1.2 and d* = 1 (corresponding 
to constant D c = 1/16). (a) Spinodal (dashed) and binodal 
(solid) lines in the (p, Do)-plane. Indicated is the behavior 
for global density p = 1: Increasing the propulsion speed Do 
(vertical arrow), the homogeneous profile loses linear stability 
when reaching the spinodal at d c . For any quenched Do, the 
coexisting densities are predicted by the points on the binodal 
(horizontal arrows for v c ). (b) The corresponding tilted free 
energy density /ad(p)—p(p—po) for Do = d c . The continuation 
of the dashed line from the upper panel indicates the inflection 
point for p = 1. 


fad (p) = /ad(Po) + P{P ~ PO ) 

+ 2 (^ >e - Yes) ^ _ Po ^ 2 + 12 < ’ 2 ^ _ P °) 4 ’ 

with 

Po = P = Po . (33) 

The bulk free energy density becomes a non-convex func¬ 
tion for speeds 

v 0 ^ 4 VD e = v *, (34) 


with the coexisting densities 


(bi) , ^ 

P± = Po ± yr 


(36) 


Consequently, the coefficient p represents the chemical 
potential, which according to Eq. (35) is equal for the two 


coexisting phases of densities p^) 1 * identifying the binodal. 
Since /" d (p) = X>(p), the inflection points defining the 
mean-field spinodal 


( s p) I 1 
P± 


Vn - Vi 


(37) 


where (p) depends on the global density. As the speed 
is increased beyond u*, the suspension enters the two- 
phase region. A common-tangent construction to mini¬ 
mize the bulk free energy density in the non-convex re¬ 
gion leads to 


d_l 

dp 



= P 


(35) 


coincide with the limit of linear stability as expected. 

The physical picture following from this discussion is 
thus that of an active suspension undergoing a phase sep¬ 
aration into a dense and a dilute phase, see Fig. [3] for an 
illustration. Although we study a system that is inces¬ 
santly driven away from thermal equilibrium, its large- 
scale evolution is apparently that of an effective equilib¬ 
rium suspension, where the speed vq plays the role of 
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an inverse temperature. In particular, the existence of a 
free energy-like functional guarantees a unique stationary 
state in which this functional becomes minimal, 


d l = [ d 2 r—^ 

d t J 6 P dt 



<0 


(38) 


inserting Eq. (291. In qualitative agreement with numeri¬ 
cal simulations "EDI ED ESI , the mean-field theory predicts 
a binodal enclosing the two-phase region and a spinodal 
within this region. The condition 


3r;* 

n = JT. =p ' 


(39) 


with = £(/!*, w*) defines the critical point (p*, v*) at 
which binodal and spinodal meet. 

To conclude this section we comment on two points: 
First, the functional Eq. (301 does not contain a gra¬ 
dient term penalizing sharp interfaces between low and 
high density phase, and we will come back to this point 
in Sec. IV B| Second, the pressure P(p) [Eq. (18)] does 
not have to be equal in the coexisting phases. The inter¬ 
face will consist of particles pointing into the dense phase 
(otherwise particles will leave the interface back into the 
dilute phase). Hence, the orientation p does not van¬ 
ish in the interface allowing for a jump of P{p) following 
Eq. (24). In addition to the pressure P(p), from a free 
energy density /(p) one can derive the “thermodynamic” 
pressure 


Pf(p) = - 


d(af) 

da 


= -f{p) + pf(p) 


(40) 


with Pf ^ P, where a is a (small) area with locally ho¬ 
mogeneous density p oc 1/a. This pressure now becomes 
equal for the coexisting densities Eq. (36). 


IV. WEAKLY NON-LINEAR ANALYSIS 


A. Dominating mode 


Employing the adiabatic solution Eq. (241 allows to 
reduce the original hydrodynamic equations to a single 
equation of motion for the density alone, which, more¬ 
over, can be cast into a form involving an effective free 
energy. We now attempt to study the behavior of the 
hydrodynamic equations in a more systematic way em¬ 
ploying a small expansion parameter [S5j . To this end we 
will consider the state points (p,v c ) along the instability 
line with force imbalance coefficients 




(41) 


following from the condition D e + a c /3 c = 0. We study 
speeds 


v 0 = v c (l +e) 


(42) 


close to the instability line with \e\ -C 1. We expand 
the coefficients a(vo) = a c + aye + 0(e 2 ) and /3{ v o) = 
(3 C + /3i£ + 0(e 2 ) into Taylor series assuming that they 
are analytic functions of the speed Vq. To leading order 
a/3 ~ —D e + aye with new coefficient 


<7i = a c f 3i + a\/3 c . 


(43) 


Hence, the fastest growing wave vector Eq. (23) behaves 
as q 0 ~ y/e and the growth rate of structures with this 
wave vector becomes a/go) ~ —aiq 2 e ~ e 2 . In the fol¬ 
lowing it will be more convenient to employ non-negative 
e ^ 0 and let ay A ±|ay| so that the sign of ay deter¬ 
mines on which side of the instability line we are: for 
ay > 0 (small) fluctuations decay while for ay < 0 the 
suspension has become unstable. 

We now aim to derive an equation of motion for the 
density fluctuations on the scale of the dominating mode. 
Since these fluctuations evolve on the length 1/go and 
grow with time scale l/a/go) we rescale length and time 
leading to 


8 t A e 2 d t , V A v'iV, 


(44) 


which we plug into Eqs. (19) and (|20|). 


B. Close to the critical point 

In a first step, we expand the local density and orien¬ 
tation as 


p = P + yfeC + £C ^ + £' 


=.3/2 c (3/2) 


p = ep (l) +£ 3/2 p (3/2) 


(45) 

(46) 


To lowest order the magnitude of density fluctuations is 
thus ~ y/s, viz. the response is Sp oc e 1 / 2 as expected 
close to a critical point with mean-field exponent |. The 
expansion form for p has been chosen to match pow¬ 
ers. Plugging all expansions back into the hydrodynamic 
equations together with Eq. ( |44| ) , we collect terms of the 
same order e. To lowest order, we find = —/3 c Vc 
and therefore 


0 = (A, + a c /3c)V 2 c, (47) 

which is fulfilled for any perturbation c(r, t) since the 
expression in brackets corresponds to the instability con¬ 
dition and thus vanishes. 

To next order we find p( 3 / 2 ) = — p c \/c^ + ( c (cVc) 
leading to 

0 = -2 gV ■ (cVc) = -gV 2 c 2 (48) 

with another coefficient 

9 = ^Cc(ay +/3 C ) = ~v 2 ^0. (49) 

For non-vanishing c this condition is fulfilled only for 
g = 0, which implies v c = v*. Hence, the expansion 
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scaled density deviation c scaled density deviation c 



scaled density deviation c 


FIG. 4: Possible shapes of the free energy density Eq. (601: (a) In the linearly unstable region for a i < 0. Due to volume 
exclusion there is a maximal c+ (vertical dashed line) at which the density has a minimum. (b,c) For eri > 0 two cases are 
possible: (b) the uniform density profile corresponding to c = 0 is metastable (the second minimum) or (c) the uniform density 
profile is globally stable. 


Eq. (45) for the density fluctuations is only valid in an 
e-environment of the critical point (/5*,u*). 

Assuming g = 0, at the next order we finally obtain an 
equation of motion 


d t c = ox V 2 c + C c V ' (c 2 Vc) - D 2 V 4 c = V 


,6Fi 


1/2 


Sc 


(50) 


for the density fluctuations. The right hand side can 
be expressed as the functional derivative of a potential 
function (an effective free energy) 


1/2 

with bulk term 


F^ /Jcl = / d 2 r 


^|Vc| 2 + / 1/2 (c) 


f 1/2 (c) = ^(TlC 2 + ^CcC 4 - 


(51) 


(52) 


In contrast to the adiabatic solution Eq. (30), the ex¬ 
pansion in £ now includes the customary square-gradient 
term. 

As the final step, we restore the density through insert¬ 
ing c = {p — p)l\fe and reverting the scaling Eq. (44). 
While the gradient term is invariant, the bulk contribu¬ 
tion to the free energy density becomes 

£ 2 /i/2 = \(*ie{p - pf + ^Cc(P - P) 4 - ( 53 ) 

In order to show that this expression for the bulk is equiv¬ 
alent to Eq. (32) obtained from the adiabatic solution, we 
first consider 


D -^l-D - —(l+£) 2 «-—£ 
e 16 6 16 1 ' E ' 8 


(54) 


up to linear order in e. On the other hand, we obtain 


o' 1 = 


(55) 


inserting the force imbalance coefficient Eq. pd] ) at the 
critical point v c = u*. This demonstrates that the adia¬ 
batic solution / a( j coincides with the result f 1 / 2 for the 


bulk free energy density of a more systematic expansion 
close to the critical point. In addition, the later route 
also yields a square-gradient term that describes the cost 
of creating density inhomogeneities. This term is impor¬ 
tant in the coarsening dynamics and allows to predict an 
interfacial tension from the mean-held theory. Note that 
such a square-gradient term can already be anticipated 
from the dispersion relation Eq. (22) from the g 4 term 
with ( a c /3 c ) 2 = D 2 . The coefficient defines a length D e 
that can be interpreted as an effective interaction range. 


C. Away from the critical point 


We found in the previous section that an expansion in 
density fluctuations of order ~ sfe only holds for v c = 
u*, i.e., at the critical point. To lift this restriction and 
to be able to move along the instability line we need 
to satisfy Eq. (48) while at the same time j / 0. To 


achieve this within the £-expansion of the density one has 
to demand that the leading contribution is of order ~ £ 
(instead of the larger ~ yfe). Effectively, this reduces the 
magnitude of density fluctuations. We, therefore, arrive 
at the expansion employed in Ref. [451 . 

p = p + ec + e 2 c^ + • • • , (56) 

p = v/ifep^ + £ 2 p (2) -I -]. (57) 

As before, we collect terms of the same order. At lowest 


order we again find Eq. (47). The next order already 


leads to the equation of motion 


,SF 


d tC = cri V 2 c - 2gV • (cVc) - D 2 V 4 c = V 2 (58) 

be 

for the density fluctuations away from the homogeneous 
profile, where 

p (2) = —/3 c Vc (2) - ft Vc + D e V 2 p (1) + CccVc (59) 

has been used. The form of the free energy functional 
coincides with Eq. (51). In particular, we obtain the 























same square-gradient term as before. However, the bulk 
free energy density now reads 


/(c) = p lC 2 - pc 3 . 


(60) 


Here the previously introduced coefficients a± and g be¬ 
come the coefficients of the square and cubic term, re¬ 
spectively. 


The striking feature of Eq. (60) is that a c 4 term sta¬ 


bilizing the dense phase is missing. Nevertheless, due to 
volume exclusion, a real suspension will not collapse into 
a single point but reach the predicted phase-separated 
state as demonstrated in experiments and simulations. 
This mechanism of volume exclusion works on length 
scales corresponding to the particle size implying c < c+ 
with c + corresponding to the maximal density. The 
damping thus arises from couplings to scales that are not 
included in the systematic expansion at the above order. 
Thus the reason for the missing c 4 term is the scale sepa¬ 
ration between the particle size and the length over which 
density fluctuations are coarse-grained. The different 
possible shapes of the free energy density /(c) depending 
on the system parameters are sketched in Fig. [4j Even 
without the c 4 term, we can derive a differential equation 
that describes the instability line, see Appendix |B| 


V. DISCUSSION 


A. Off-critical quenches 


The two local free energy densities / 172 (c) and /(c) 
for the two different scalings can be derived from the 
same global free energy. To understand this and the 
role of the asymmetric c 3 term in the free energy den¬ 
sity Eq. (601 it is instructive to recall the situation (for 
passive suspensions) with a symmetric free energy den¬ 
sity f(ip) = |v ? 2 + jif 4 with order parameter <p oc p — p*. 
The Cahn-Hilliard equation models the evolution after 
a quench from a stable homogeneous density ip = — Cq 
into the two-phase region, where c$ = 0 implies a quench 
through the critical point. For an off-critical quench, we 
set ip = c — Co with c the deviations away from the ini¬ 
tial homogenous density as before. The Cahn-Hilliard 
equation then reads 

d t c = V 2 [(a + 3 kCq)c — 3 nc 0 c 2 + nc 3 ] — D 2 V 4 c (61) 

dropping constant terms which vanish due to the spa¬ 
tial derivative. Cleary, this result becomes Eq. (58) with 
g = 3kco and o\ = a + g 2 /(3n) when adding a repul¬ 
sive term \kc 4 to /(c). The coefficient a(v 0 ) should be a 
monotonically decreasing function of speed Vq . It changes 
sign at vq = u* and becomes negative for vq > v *. A given 
global density p determines v c = v c (p) and g = g(v c ) 
[Eq. (49)]. Exactly on the spinodal oq = 0 holds, which 
implies a = —g 2 /{3 k) for the quadratic coefficient of the 
global free energy. Note that at v c = v * we have g = 0 


and, therefore, we recover the function / 172 (c) [Eq. (52)] 
if in addition we set n(v*) = Hence, /(c) extends 
the lower order solution / 172 (c) to higher speeds Vq but 
with undetermined coefficient k( 7 j 0 ), which does not fol¬ 
low from the e-expansion. 


B. Nucleation behavior close to the spinodal 


The Cahn-Hilliard equation has been derived origi¬ 
nally to describe spinodal decomposition , the homoge¬ 
neous, barrierless onset of phase separation throughout 
the system in response to a quench beyond the spinodal 
bounding the instable region. Interestingly, for the spin¬ 
odal itself the Cahn-Hilliard equation predicts a change 
from a continuous to a discontinuous transition (see, e.g. 
Ref. S3). 


To discuss this effect in the present context of active 
Brownian particles, we symmetrize and rescale the free 
energy density Eq. (60) (including the repulsive c 4 term) 
leading to the scaling functions 


/±fa)= r ^[/(c)-/(c 0 )] 

m\ 

= p±v + ^(±1 - i> 2 + p 4 


(62) 


with scaled density = \J k/|cti |(c — Co), where Co = 
g/(3ft) as in the previous section. These functions depend 
on the single parameter 


r = 


gc 0 

M 


9 

3k|cti| 


0 


(63) 


combining the three coefficients a ±, g, and k, with irrele¬ 
vant p.±(T) akin to a chemical potential. Corresponding 
to the sign of 01 , /_ describes the effective free energy 
density in the unstable and /+ in the linearly stable re¬ 
gion. 

The parameter T is defined along the spinodal. In par¬ 
ticular, the critical point (p*,v*) corresponds to T = 0 
and going away from this point along the spinodal, the 
value of T increases. In Fig. [5] the resulting “phase dia¬ 
gram” is shown in the ( 77 , T) plane. Due to the scaling 
employed in the previous section, the functions f±{ij) 
correspond to fixed (small) e and we now need to distin¬ 
guish the two cases shown in Fig. [5](a) and (b) depending 
on which side of the spinodal the system resides. Shown 
are the binodals calculated from the common tangent 
construction and the spinodals corresponding to the in¬ 
flection points, cf. Sec. |III B| The uniform density corre¬ 
sponds to setting c = 0, i.e., g = —y/T/3. For o\ > 0 the 
uniform density is linearly stable but becomes metastable 
for T > | (it crosses the binodal). This implies that for 
r < | crossing the spinodal (jumping from <7\ > 0 to 
(7 1 < 0) corresponds to a continuous transition. In con¬ 
trast, for r > | the transition becomes discontinuous. 

It is instructive to also consider the bifurcation dia¬ 
grams shown in Fig. [5](c,d). Here we show schematically 
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cq > 0 (linearly stable) cq < 0 (unstable) 



scaled density r| scaled density r| 



FIG. 5: Schematic “phase diagrams” obtained from the scal¬ 
ing functions (a) /+ and (b) /_, see Eq. ( |62|) . Shown are co¬ 
existing densities (binodals, outer solid lines) and the limits 
of linear stability (spinodals, inner dashed lines). For eri > 0 
the uniform density (black line) is linearly stable (it becomes 
metastable for F > |) and for o\ < 0 it is always unstable. 
(c,d) Sketch of the bifurcation diagrams showing the devia¬ 
tion from the uniform density (the “amplitude”, solid lines are 
linearly stable and dashed lines unstable) vs. the control pa¬ 
rameter e for fixed T. The diagrams of (a) and (b) correspond 
to ±|e| on both sides of the bifurcation point as indicated by 
the connected dots for two representative values of V: (c) Su¬ 
percritical pitchfork bifurcation for F < | corresponding to a 
continuous transition, (d) Subcritical bifurcation for T > | 
corresponding to a discontinuous transition. 


the solutions of the amplitude equation, which corre¬ 
spond to the extrema of the effective free energy. Plotted 
is the density with respect to the homogeneous density 
as a function of e = vq/v c — 1 for fixed P. For any |e| > 0 
close to the instability at e = 0, the (scaled) solutions 
can be directly read off Fig. [5](a) for e < 0 (ex, > 0) and 
Fig. |5](b) for e > 0 (a± < 0) as indicated for two values of 
r. These correspond to a continuous and discontinuous 
transition, respectively. 


C. Non-local speed 

Within the weakly non-linear analysis an integrable 
square-gradient term stabilizing domains appears nat¬ 
urally. To study phase separation kinetics, Cates and 
coworkers have followed a different route and posit 
that the active Brownian particles sample the density 
on a length scale A larger than the interparticle spac¬ 
ing EE :2] ■ To lowest order the speed then becomes a 
non-local function v(p) with p = p + A 2 V 2 p such that 

v(p) « v(p) + v'(p)X 2 X/ 2 p. (64) 


Plugging such a non-local speed into the adiabatic solu¬ 
tion Eq. (241 produces the desired square-gradient term 


but, strikingly, would lead to an equation of motion for 
the density that also involves non-integrable terms, i.e., 
it is not longer representable as the functional derivative 
of an effective free energy. Apparent consequences are 
discussed in Ref. El 


Replacing the speed Eq. (13) with the expression 
Eq. (64), we can again study the systematic e-expansion 
of the hydrodynamic equations (15) and close to the 
instability line. While the lowest order of the e-expansion 
remains unchanged, the non-local speed modifies the ex¬ 
pression Eq. (59) 


d 2 ) 


J2) 


1 


^ + 2^ cA V ( Vc ) 


(65) 


for the orientation, where we have used v'(p) = —f. Con¬ 
sequently, the bulk free energy density Eq. (601 remains 


the same and the only effect is that the coefficient of the 
square-gradient term is replaced by 


D 2 i->- D 2 + -a c p( c X 2 . 


( 66 ) 


Since the additional term is positive, this corresponds to 
a larger interaction range as one would intuitively expect. 

In summary, non-integrable terms were introduced 
when employing the adiabatic approximation. They do 
not appear in a systematic treatment. Instead, the sys¬ 
tematic treatment naturally produces the well-known 
squared density gradient term in the effective free energy 
functional. As a consequence, close to the instability line 
and at onset (mean-field) phase separation kinetics of 
active suspensions is predicted to not qualitatively dif¬ 
fer from that of passive suspensions even for a non-local 
speed. 


D. Numerical simulations 

Although there are already quite a few numerical stud¬ 
ies of two-dimensional active Brownian particles |EB- 
33. 43, 51 , the determination of the full phase diagram 
is still an open task. As shown, the mapping to a free en¬ 
ergy strictly holds only close to the instability line. More¬ 
over, mean-field treatments are known to yield incorrect 
quantitative predictions close to critical points due to a 
diverging correlation length and the ensuing large fluc¬ 
tuations. These fluctuations make numerical sampling 
notoriously difficult in the vicinity of a continuous tran¬ 
sition. This problem is already present in passive sus¬ 
pensions but appears to be even more severe in active 
suspensions due to two reasons: (i) fluctuations are even 
larger and (ii) advanced sampling techniques for systems 
out of equilibrium are not (yet) available. 

We examine the data from Refs. Jij, (43] 50] ob¬ 
tained from particle-resolved simulations of discoid self- 
propelled disks with diameter a interacting via the purely 
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FIG. 6 : (a) Coexisting densities for a model suspension, see text for details. Shown are (•) the values from Bialke et al. [50] 
and (□) Redner et al. !JT]. The dashed lines are fits of Eq. (681 with (o) indicating the extrapolated position of the critical 


point, (b) Fitted binodal and critical point as in (a) together with the average fraction m of particles in the largest domain 
(scale bar on the right, from Ref. 32) • Only state points with m ^ 0.1 are shown, (c) Comparison of m at if> = 0.5 for two 
methods: (o) from a cluster analysis to identify dense domains and (•) from the coexisting densities via Eq. (691. 


repulsive Weeks-Chandler-Andersen potential [52] 

B(r) J*[(«W"- ( »W‘] + < (r<*£0 ( 67 ) 

W [0 (r > 2 1 / 6 cr) V ' 

with potential strength e. We set a = 2~ 1 ^{a/£) so 
that repulsive interactions are only present for overlap¬ 
ping particles. In the following, we will present numeri¬ 
cal results using the global area fraction </> = ( a/£) 2 np/A , 
which is simply proportional to densities. In principle, 
the rotational diffusion coefficient D r is a free parameter 
of the model. However, for self-propelled colloidal parti¬ 
cles it is reasonable to assume (see, e.g., supplementary 
materials of Ref. H that the no-slip boundary condition 
still holds, which implies D r = 3Do/a 2 . This additional 
assumption relates the two lengths a = V3£. For bet¬ 
ter comparison, in the following we report the original 
numerical values for the dimensionless speed Vq = v / 3r , o- 

The accurate determination of coexisting densities in 
numerical simulations is subtle. Redner et al. have ex¬ 
tracted these densities from histograms [31] . In Ref. ED] 
a different strategy is followed, in which a non-square 
simulation box with periodic boundaries is employed. In 
such a geometry a “slab” of the dense phase forms, which 
aligns with the shorter box edge. This allows a very accu¬ 
rate determination of density profiles and the coexisting 
densities. The result for N = 4,900 disks at global area 
fraction <fi = 0.5 using e = lOOfeeF is shown in Fig. |6][a). 
We have obtained data for N = 10, 000, which agree with 
the coexisting densities for the smaller system demon¬ 
strating that there is no finite-size dependence. Also 
shown is the data from Ref. SB who have used e = ksT 
and N = 15,000. The low-density branches agree very 
well while for the high densities the interaction strength 
plays a more important role. Clearly, the system with 
the higher e is less compressible. Still, both studies agree 
qualitatively. Note that below v 0 < 60 fluctuations are 
so strong that a reliable determination of coexisting den¬ 


sities is not possible anymore. 

We have fitted the coexisting densities (the binodal) of 
Ref. [50] with the functional form 


0 * 


1 + a± 



+ b ± 



( 68 ) 


inspired by the famous Guggenheim plot [53]. Here, a± 
and b± are fit parameters, and (f>± are the coexisting area 
fractions. This fit works surprisingly well over a wide 
range of speeds and allows to extrapolate the position 
of the putative critical point: ~ 0.62 and D* — 48. 

In Fig. [ 6 jb) , the fitted coexisting densities are overlaid 
by results from instantaneous quenching: The passive 
suspension (vq = 0 ) is equilibrated at a given global den¬ 
sity. The suspension is then quenched to the final speed 
Vo and relaxed to the steady state. A cluster analysis 
is performed for the recorded configurations to identify 
dense domains of N c particles, from which the averaged 
fraction m = ( N c )/N of particles in the largest cluster 
( i.e ., dense domain) is extracted. Note that only state 
points with m ^ 0.1 are shown, i.e., the largest domain 
contains at least 10 % of all particles. 

Fig. [ 6 ](b) is compatible with the proposed mapping to 
passive phase separation. Below area fraction cf> ~ 0.3 
no spontaneous phase separation is observed, which thus 
roughly indicates the location of the spinodal. There is 
a larger region around the putative critical point where 
spontaneous phase separation seems to occur even below 
the extrapolated binodal, which, however, could be in¬ 
terpreted as a finite-size effect. In a finite system, large 
fluctuations appear as dense domains but would not lead 
to full phase separation in a larger system. However, this 
does not explain the observed phase separation at higher 
area fraction <f> ^ 0.7. At these high densities details of 
the interparticle interactions and the associated particle 
length scale may become important, which is not cap¬ 
tured by the coarse-grained point of view of the weakly 
nonlinear analysis. 













11 


Previously in Ref. [52] we have reported data that 
would suggest a continuous transition (for cf> = 0.4 and 
(f> = 0.5). In line with this observation, in Ref. [43] an 
abrupt onset of hysteresis for area fraction cf) < 0.32 has 
been observed, which thus agrees with the mean-field pre- 
dicti on of a change from continuous to discontinuous (see 
Sec. VB|). Whether this observed change is a numerical 


artefact due to the vicinity to a single critical point (as 
in passive phase separation) or a genuine non-equilibrium 
effect remains to be investigated. 

For completeness, we note that for a single dense do¬ 
main (as observed at intermediate system sizes) a simple 
relation between m and the coexisting densities exists via 
the lever rule, 


m{v 0 ) 


1/0- - V0 

1 /0- - 1/0+ ' 


(69) 


This is demonstrated in Fig. [Gjc) . The result for m 
obtained from the cluster analysis is slightly larger, 
which may be attributed to the interfacial particles being 
counted towards the dense phase. 


VI. CONCLUSIONS 

To summarize, we have followed a systematic route 
from the microscopic dynamics to the large-scale Cahn- 
Hilliard equation to address the phase behavior of ac¬ 
tive Brownian particles. In Ref. [32], starting from the 
microscopic dynamics, effective hydrodynamic equations 
have been derived for the temporal evolution of (weakly 
perturbed) density and average orientation. Two param¬ 
eters enter these equations: the free speed vq and the 
global density p. All details of the particle interactions 
are contained in the single function £ = f(p,v o), which 
quantifies the force imbalance due to the interplay be¬ 
tween repulsive forces and directed motion. 

The hydrodynamic equations exhibit a dynamical in¬ 
stability at v c = v c (p), which can be determined through 
a linear stability analysis [301 EH GSj . To go beyond the 
linear regime, we have performed a weakly non-linear 
analysis [44] close to the instability line using as the small 
expansion parameter e = vq/v c — 1. This approach yields 
an equation of motion for the density alone, which is 
formally equivalent to a Cahn-Hilliard equation. It in¬ 
volves a local effective free energy. We have discussed 
two expansions holding for different ranges of the speed 
v c , which together agree with the scenario of mapping 
active to passive phase separation close to the loss of 
linear stability. What happens further away from the in¬ 
stability line? Including higher orders of the expansion 
is of course possible but will lead to a large number of 
additional terms. Already at the next order the non¬ 
equilibrium nature of the active suspension will become 
manifest since a description in terms of density fluctua¬ 
tions alone is not possible anymore and the time evolu¬ 
tion of the polarization has to be included. Hence, there 
seems to be no real advantage in going to higher orders 


since one can also study the original hydrodynamic equa¬ 
tions. 

One should stress that we have described a mean- 
field scenario for the one-point density. Stationary two- 
point correlations are an input to the theory and higher 
order correlations are neglected. Gaussian noise could 
be added, which accounts for uncorrelated fluctuations 
but of course does not restore the missing correlations. 
For passive suspensions it is well known that mean-field 
free energies have to be regularized by the Maxwell con¬ 
struction to be thermodynamically valid. Moreover, ex¬ 
perimentally and in simulations no sharp loss of linear 
stability is observed and the “spinodal” line is a pure 
mean-field concept (54] . Still, one might wonder why 
for active Brownian particles a mean-field description 
nevertheless seems to reproduce domain morphologies of 
particle-based simulations to such a good degree [39]. It 
thus remains to be tested numerically to which extend 
the theory outlined here is valid. Some evidence has been 
discussed in Sec. IV Dl but more detailed numerical inves¬ 
tigations are clearly needed. 

Despite recent progress there are many open questions 
even for the simple minimal model studies here. For ex¬ 
ample, the determination of the full phase diagram is 
still an open issue both theoretical and numerical. As we 
have shown here, in the mean-field picture already the 
instability line has a richer structure than anticipated 
previously. It will be particularly interesting to study 
speeds close to u* in order to determine the nature of the 
disorder-order transition. However, simulations in this 
parameter region are hampered by the large critical fluc¬ 
tuations and will require a more detailed study of finite- 
size effects than available so far. Another open issue is 
the relation between effective free energy and mechanical 
pressure [SOI EH ESHSZI • 

Finally, while here we have studied analytically the 
simplest version of active Brownian particles (namely 
monodisperse repulsive disks) in two dimensions with¬ 
out alignment and without hydrodynamic interactions, 
we note that multiple extensions of this model have been 
discussed in the literature: mixtures of active and passive 
particles [53], role of attractive interactions [53 EU] and 
alignment [61| , and polydispersity in connection with the 
glass transition [521 <13 ■ 
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Appendix A: Run-and-tumble particles 


The equivalence of active Brownian particles and run- 
and-tumble particles on the level of the hydrodynamic 
equations has been discussed in Ref. [37. For complete¬ 
ness, we briefly sketch the derivation for run-and-tumble 
particles following the route taken in Sec. [TT] Again, we 
consider N identical particles moving in two dimensions. 
A particle moves with constant speed vq along its orien¬ 
tation e = (cos <p, sin ip) T (the “run”). After an expo¬ 
nentially distributed random run time with mean r r , the 
particle “tumbles”, i. e. it picks a random new orienta¬ 
tion. The evolution of the one-point density -0i (r, ip,t) is 
thus given by 

dt'tpi = -V • [n 0 F + voeipi] - ~ipi + -p , (Al) 


where F is the force due to interactions with other par¬ 
ticles, no the bare mobility, and p( r, t) the local num¬ 
ber density [Eq. (15)]. The second term describes the 
“death” of particles with a given orientation ip and the 
third term their rebirth with uniformly distributed orien¬ 
tation. Dimensionless quantities are introduced through 
measuring energies in units of e, rescaling time t H > T r t , 
and length r^><r with t = y/epor^. 

As a closure, we decompose the force F w — pC e '0i — 
_D e V/’i, where the first term captures the force imbalance 
along the orientation slowing down the particles. Since 
this force is not exactly aligned with the orientation there 
will also be an “evasive” motion, leading to an effective 
diffusion on larger scales. This is described through the 
second term. Inserting the force into Eq. ( |A1| ), we thus 
arrive at [cf. Eq. @] 


d t ipi = -V ■ [v(p)e - D e V]ipi -i/j 1 + 77 - (A2) 

Z7T 


with effective speed v(p) = vq — pC- It is now straight¬ 
forward to derive the hydrodynamic equations (15) and 
(17), which demonstrates the equivalence of active Brow¬ 
nian particles and run-and-tumble particles (at least for 
weakly perturbed orientations). 


these cannot be completely independent. Consider a 
global density pi. Quenching the system to a speed 
Vq = t> c (pi)(l + £ ) past the limit of linear stability, sepa¬ 
ration into dense and dilute regions will occur. The range 
of unstable scaled density deviations is cf! P ^ < c < c+ p \ 
where c^ p ^ = <ti/(2 g) < 0 is given by the inflection point 
of /(c) (the maximal value is irrelevant for the fol¬ 
lowing argument). Clearly, a suspension with global den¬ 
sity P 2 = Pi + ec^! P ' 1 should then become unstable beyond 
the speed vo = v c (p 2 ). We can exploit this consistency 
condition to derive a differential equation that describes 
the instability line, see Fig. [7] for a sketch of the deriva¬ 
tion. 

To this end, consider the derivative 


dp 

dv 


Vl 


y P 2 - Pi 
lim - 

V2-M>1 V2 — Vl 


lim 

v 2 —MJi 


£C 


(sp) 


V 2 - Vl 


(Bl) 


at speed iq = v c (p±). Inserting e = ( V 2 — vi)/vi, we 
obtain 


_ oj_ 

dv 2gv 


(B2) 


Knowing one point on the instability line, we can inte¬ 
grate this differential equation and invert the solution to 
obtain the full instability line v c (p), see Ref. [43] for an 
example. 



Appendix B: Differential equation for instability line 


The effective free energy density Eq. (60) is a local 


function that describes the scaled density deviations in 
the vicinity of a point on the instability line. Two points 
along this line have different coefficients <ti and g. Still, 


FIG. 7: Sketch for the derivation of Eq. (B21 for the instability 


line (solid line). F or a global density pi the bulk free energy 
density /(c) [Eq. (60)] is sketched in the inset, whereby the 
homogeneous density corresponds to c = 0. Densities down 
to P 2 = pi + ec)! P ' ) are unstable, which thus corresponds to a 
second point on the instability line. 
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